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This paper uses the assumptions of ergodicity and a microcanonical distribution to compute es- 
timates of the largest Lyapunov exponents in lower-dimensional Hamiltonian systems. That the 
resulting estimates are in reasonable agreement with the actual values computed numerically cor- 
roborates the intuition that chaos in such systems can be understood as arising generically from a 
parametric instability and that this instability can be modeled by a stochastic-oscillator equation 
(c/. Casetti, Clementi, and Pettini, Phys. Rev. E 54, 5969 (1996)), linearised perturbations of a 
chaotic orbit satisfying a harmonic-oscillator equation with a randomly varying frequency. 



PACS number(s): 05.45. -l-h, 02.40.-k, 05.20.-y 
I. INTRODUCTION AND MOTIVATION 

By definition, Lyapunov exponents probe the aver- 
age linear instability of chaotic orbits in an asymptotic 
t ^ GO limit . Their computation thus involves solving 
a matrix harmonic-oscillator equation with characteristic 
frequencies that vary in time. In the context of a geo- 
metric description - which is convenient but by no means 
necessary -, this equation can be reinterpreted 
cobi equation [i.e., equation of geodesic deviation) for 
motion in an appropriately defined curved space, e.g., by 
introducing the Eisenhart metric ||^ . 

It has been long known |^,^ that geodesic flows in a 
space with everywhere negative curvature are unstable in 
the sense that nearby orbits diverge exponentially; and, 
for this reason, there was an implicit assumption in much 
earlier work that chaos could often be understood as a 
manifestation of negative curvature. However, as em- 
phasized by Pettini [|j , in many systems chaos cannot be 
attributed to negative curvature. In many cases, the av- 
erage curvature is positive; and indeed, there are many 
known examples of nonintegrable Hamiltonian systems 
{e.g., the finite order truncations of the Toda |^ poten- 
tial) which admit large measures of chaos even though 
the curvature is everywhere positive. The curvature as- 
sociated with the Eisenhart metric can be negative only 
if the second derivative of the potential becomes nega- 
tive. Instead, it would seem natural to understand chaos 
as refiecting a parametric instability. 

The Jacobi equation for a regular periodic orbit re- 
duces to a multi-dimensional Hill equation, i.e., a 
harmonic-oscillator equation with frequencies that ex- 
hibit a periodic time dependence. For certain amplitudes 
and periodicities, the solutions to such an equation re- 
main bounded (or grow at most as a power law in time), 



this corresponding to stable periodic orbits. However, 
for other amplitudes and periodicities, the solutions grow 
exponentially, this corresponding instead to an unstable 
periodic orbit 

Since chaotic orbits are aperiodic and (in some sense) 
'random,' one might instead suppose that one can model 
the Jacobi equation describing a linearised perturbation 
of a chaotic orbit as a stochastic harmonic-oscillator equa- 
tion, in which the time-dependent frequencies vary in a 
random fashion. Given this assumption, the key issue be- 
comes one of identifying the stochastic process, i.e. the 
form of the colored noise, which can capture correctly 
solutions to the Jacobi equation. 

If, for fixed potential and energy, almost all of the con- 
stant energy hypersurface is chaotic, as is true generi- 
cally for D > 2 (provided that the energy E is the only 
time- independent constant of the motion), it would seem 
reasonable to infer that the parameters for the oscilla- 
tor equation should be estimatable assuming ergodicity. 
What this means is that one can assume an invariant 
measure corresponding to a uniform population of the 
constant energy hypersurface, i.e., a microcanonical dis- 
tribution. If, furthermore, one is concerned with compar- 
atively high-dimensional systems, the computationally 
awkward description in terms of a microcanonical distri- 
bution can be replaced by a more user-friendly descrip- 
tion based on a canonical distribution: In the spirit of or- 
dinary thermodynamics, one can argue that the canonical 
and microcanonical ensembles should yield nearly iden- 
tical results in the large D limit. 

Given this logic, Casetti, Clementi, and Pettini Q de- 
veloped a 'thermodynamic' theory of chaos which they 
used to obtain very good estimates of the values of the 
largest Lyapunov exponents for two well studied phys- 
ical systems. To do this, they: i) extracted from the 



1 



full _D-dimensional Jacobi equation an 'isotropized' one- 
dimensional oscillator equation which they argued should 
capture the chaotic behavior of typical orbits; (ii) derived 
the statistics of their assumed stochastic process in the 
context of a canonical ensemble description; and then 
(iii) showed that, for two seemingly generic models, solu- 
tions to the resulting equation yield reasonable estimates 
of the largest Lyapunov exponent, at least for D > 100 
or so. 

An obvious question is whether this logic can also be 
exploited to provide reasonable estimates of the largest 
Lyapunov exponent for lower dimensional systems, say 
D = 2 or D = 3. As discussed in the concluding Sec- 
tion, there are a variety of settings where it would be 
convenient if one could estimate these values without re- 
sorting to detailed numerical computations. Arguably, 
however, this is not the most important point. Rather, 
the foremost objective is to implement a simple physi- 
cal picture of the origins of chaos in lower-dimensional 
Hamiltonian systems. To the extent that the Casetti et 
al proposal, or some variant thereof, can provide reason- 
able estimates of Lyapunov exponents in these systems, 
one would seem justified in visualising chaos as arising 
from a parametric instability manifested by a stochastic- 
oscillator equation. In other words, one will have a clear 
new paradigm in terms of which to interpret the origins 
of chaos in lower-dimensional Hamiltonian systems. 

II. AN ILLUSTRATIVE EXAMPLE 

The validity of the formula for the largest Lyapunov ex- 
ponent derived by Casetti et al. was tested for one simple 
toy model. The model is motivated by recent observa- 
tions of elliptical galaxies, which suggest that these ob- 
jects can exhibit significant deviations from axisymmetry 
and that they often have a high-density cusp at their cen- 
troids, perhaps associated with the presence of a super- 
massive black hole. The stars in a real galaxy populate a 
6Af-dimensional phase space, with N denoting the num- 
ber of stars in a system. Considering that fine structure 
due to localised irregularities and granularity will take 
a long time to manifest itself, it is of interest to model 
the system in terms of its coarse-grained 6-dimensional 
phase space in expectation that the time scale associ- 
ated with the coarse-grained potential will represent the 
shortest time scale for macroscopic evolution. Thus, a 
model potential for such systems comprises the sum of 
an anisotropic harmonic potential and a spherical Plum- 
mer potential: 

V{x, y, z) = \{a^x^ + + c'z') - (2.1) 

withr^ = x'^+y'^ + z'^, = 1-A, \p- = 1, and = 1 + A. 
A parameterizes the ellipsoidal geometry, and e functions 
as a "softening parameter" which is set at e = 10~^ for 
numerical simulations. 



The theory of Casetti et al., described in Section IV 
A below is analytic, and within this formalism e acts as 
a "free parameter" that reflects uncertainty about the 
detailed dynamical properties of the phase space. One 
knows a priori that far from Mbh the potential it is ap- 
proximately quadratic in the coordinates, and close to 
Mbh it is approximately spherically symmetric; the or- 
bits are accordingly quasiregular in these regions wherein 
there will be almost no chaotic mixing. The theory 
correctly predicts zero chaotic mixing in a harmonic- 
oscillator potential, thereby incorporating the former 
circumstance, but it also incorrectly predicts nonzero 
chaotic mixing in the spherically symmetric Plummer po- 
tential that dominates near Mbh- Thus, a nonzero e 
"regularizes" orbits near the black hole. In view of these 
considerations, the value of e used in the theory was cho- 
sen by requiring the magnitude of the harmonic poten- 
tial to be comparable to a tenth of that of the Plummer 
potential at distances "r — e" from the centroid. The 
specific choice is e = 0.5M^^^. 

FIGURE g compares the "true" Lyapunov exponents, 
computed via numerical simulations with estimates of the 
largest Lyapunov exponent derived using the Cassetti et 
al. formalism, Eqs. (4.12), (4.15), and (4.16) below. The 
numerical studies are described in detail in Ref. and 
the numerically generated curves derive from FIG. 5 in 
that paper. The figure shows how the Lyapunov expo- 
nents pertaining to chaotic orbits scale against black-hole 
mass and total particle energy E. Interestingly, the ana- 
lytic results agree closely with the numerical results, par- 
ticularly for intermediate-to-small values of Mbh- The 
agreement is still reasonable for large values of Al bh (val- 
ues that are in fact unphysically large) , though the degree 
of agreement is less good. This is as expected in that 
a black-hole mass that is comparable to the ellipsoidal 
mass will establish sizeable regions of regularity over the 
constant-energy hypersurface, and the fraction of chaotic 
orbits will be correspondingly lower . 

FIGURE ^ compares, for fixed Mbh = 0.1, the numer- 
ical and analytic Lyapunov exponents versus cUipticity 
as parameterized by A. Again, the analytic technique is 
seen to yield reasonable estimates provided A is not too 
small. As A decreases to zero, the potential approaches 
spherical symmetry and is thereby integrable, support- 
ing only regular orbits. Inasmuch as the fundamental as- 
sumption underlying the Casetti et al. formalism is that 
a substantial fraction of the orbits is globally chaotic, 
the formalism clearly breaks down for spherical symme- 
try. As discussed in Ref. |^, the fact that the numerical 
curves exhibit a great deal of structure not manifested by 
the analytic predictions reflects the fact that the phase 
space associated with the potential (2.1) is dominated by 
resonances with frequencies a, 6, and c associated with 
the harmonic contribution which are completely indepen- 
dent of initial conditions. 

The results of FIGS. Q and ^ suggest that the 6- 
dimensional phase space governed by the toy poten- 



2 



tial (2.1) exhibits global chaos and associated rapid ir- 
reversible mixing over the bulk of the parameter space. 
Can the same be said for a lower-dimensional analog, 
i.e. one corresponding to the toy potential in which 
^ = Pz = 0? FIGURE H, which provides the same in- 
formation as FIG. ^ but now for a 4-dimensional phase 
space, hints at the answer. One now sees the agree- 
ment between the numerical and analytic results to be 
less good, as would be expected because the fraction of 
globally chaotic orbits is generally much reduced over 
the 6-dimensional case. Nonetheless, the results are still 
comparable within a factor of two. 



III. THE SCOPE OF THIS PAPER 

The obvious question is whether the striking agreement 
between theory and numerics described in the preceding 
section is simply a fortuitous accident, or whether it is 
in fact generic. Can the Casetti et al. analysis provide 
reasonable estimates of the largest Lyapunov exponent 
for generic lower-dimensional Hamiltonian systems? 

Related to this is another important question: To what 
extent are the assumptions implemented by Casetti et al 
justified for lower-dimensional systems? To the extent 
that they are not justified, one might expect either (i) 
that the final formula for x which they derived is compar- 
atively insensitive to (some of) the assumptions and/or 
(ii) that modifying these assumptions might lead to im- 
proved estimates. 

These questions were addressed by a detailed explo- 
ration of orbits in the potentials discussed in Section II, 
as well as three other (classes of) potentials: 
1. The sixth order truncation of the Toda lattice a 
familiar two-dimensional potential: 



v{x, y) = ^ {x^ + y^) + x^y - 



+ x^y + ^x^y^ - + l^-e + ^iy2 ^ ^^2yi ^ 11^6^ 

(3.1) 

2. A multi-dimensional generalisation of the dihedral po- 
tential | [lO| |, for one particular set of parameter values, 
allowing for _D = 2 through D ~ 6: 



D / D \^ D 

y(,r,..,.z.) = -E^f + i E^n -I E 

1=1 \i=i J i<j=i 



2 2 



(3.2) 



3. A generalisation of the Fermi-Pasta-Ulam {FPU) (3 
model with 



V{qi, ..,qD) = E 



^ {qi+1 - qdf + (gj+i - qif 



with (7_D+i = allowing for £) = 3 through D = 6 (the 
special case = 2 is integrable). For a = 1, (3.3) reduces 
to the standard FPU model; for a < 0, the potential 
admits extrema that are local maxima, so that the local 
mean curvature can become negative. The case a = 1 
and & = 0.1 was considered by Casetti et al. for much 
larger values of D. 

Section IV of this paper begins by providing a terse 
mathematical summary of the formalism introduced by 
Casetti et al. to estimate the values of the largest Lya- 
punov exponent in higher-dimensional systems. This 
mathematical structure is then restated in much simpler 
physical language and the resulting reformulation is used 
to suggest how their analysis could be reformulated for 
lower-dimensional systems. Section V summarises the re- 
sults of extensive simulations in the potentials (2.1) and 
(3.1) - (3.3) which were used to test the validity of the 
original assumptions. Section VI then turns to the ac- 
tual values of Lyapunov exponents estimated using this 
general approach, considering both the 'true' Lyapunov 
exponents, defined in at oo limit, and short-time Lya- 
punov exponents appropriate for orbit segments of 
comparatively short duration. Estimates of the latter 
for a variety of different orbit segments evolved in the 
same potential with the same energy reveals an impor- 
tant point: Even when the estimated short-time expo- 
nents Xest differs from the 'true' exponents Xnum com- 
puted numerically by as much as a factor of two, their 
values tend to be strongly correlated. For example, orbit 
segments for which Xest is especially small correspond in 
general to orbits for which the numerical Xnum is also 
especially small. In this sense, it is clear that, even if 
the Casetti et al. formula for x is not completely sat- 
isfactory, it does capture some important aspects of the 
flow. Section VII concludes by summarising the princi- 
pal conclusions and discussing potential implications and 
extensions. 



IV. CHAOTIC MOTION AS A STOCHASTIC 
PARAMETRIC INSTABILITY 

A. The proposal of Casetti, Clementi, and Pettini 

The starting point is the reformulation of a time- 
independent Hamiltonian system as a geodesic flow in 
an appropriately defined curved space. This can be done 
in a variety of different ways, the best known of which 
involves implementing Maupertuis' Principle p^ , which 
leads to the Jacobi metric. However, from a practical 
perspective, the most convenient choice is to work with 
the Eisenhart metric 

Given a D degree-of-freedom Hamiltonian system char- 
acterised by a Lagrangian 
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aijq'q^ 



V{q\. 



..,q 



(4.1) 



(3.3) 
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with motion defined on some manifold M, consider 
the extended manifold M x R^, with coordinates 



, q , g + ), and introduce the Eisenhart metric 



= aijdq'dq^ - 2V{q)dq°dq° + 2dq°dq°+\ 



(4.2) 



Setting 9° = t and q°+^ = t/2 - dt'L(q, q) yields 
ds^ = dt^. Without loss of generality one can set 
(i-ij = Sij, making the kinetic energy a sum of quadratic 
contributions (g*)^/2, in which case the geodesic equa- 
tions reduce to Newton's equations of motion. Corre- 
spondingly the Riemann tensor simplifies greatly; its only 
nonvanishing components are 



dq^dq^ 



(4.3) 



and the Jacobi equation for a linearised perturbation be- 
comes 



0, 



(* 



1, 



(4.4) 



Were the Riemann components entering into Eq. (4.4) 
were everywhere negative, an arbitrary perturbation 
would always grow exponentially fast. Everywhere nega- 
tive curvature implies chaotic behavior and positive Lya- 
punov exponents The important point, however, 

is that, because of parametric instability, one can have 
chaotic orbits with positive Lyapunov exponents even if 
the curvature is everywhere positive. If, following Casetti 
et ai, one assumes that the curvature varies 'randomly' 
along a chaotic orbit, eq. (4.5) reduces to a stochastic- 
oscillator equation of the form 



(4.5) 



where the matrix k^j is characterised completely by its 
statistical properties. However, it is well known that, 
even if k''^ is positive definite for all times, can grow 
exponentially fl^ . 

Especially in high dimensions, matrix equations be- 
come difficult to solve either numerically or analytically. 
For this reason, Casetti et al. proceeded by replacing this 
D-dimensional equation with a simpler one-dimensional 
equation which aims to capture its 'average' behavior. 
Formally, they start by observing that the Riemann ten- 
sor can be decomposed into two pieces: 
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RjkSlj 



W 



jki 



(4.6) 



where Rij the Ricci tensor and VFjfe; is the Weyl pro- 
jective tensor. For the case of an isotropic space, the 
Weyl tensor vanishes identically and Rjicf q^ / {D — 1) re- 
duces to the (constant) sectional curvature. The crucial 
assumption then is that, even though the space is not 
isotropic, it should appear nearly isotropic when viewed 



on comparatively large scales. In the context of such a 
'quasi-isotropic' approximation, R^^^ and I fixed this: 
ui'^ = ^^''^^jik become diagonal, and Eq (4.6) reduces to 
a one-dimensional equation of the form 



d^ 

df^ 



k{t)C = 0. 



(4.7) 



All that remains is to specify the statistical properties of 
the random process k{t). 

For a generic Hamiltonian system, the form of k could 
be quite complex. However, given the assumption that 
D is large, one might expect that the complicated details 
will largely wash out. Thus, Casetti et al. argue in the 
spirit of the Central Limits Theorem that the curvature 
fiuctuations in different directions can be approximated 
as nearly independent and, at any instant, Gaussianly 
distributed. It then follows that 



d^ 
df^ 



(4.8) 



where, in terms of the quantity K = d^V / dq^dqi = AT^, 
which has mean {K) and dispersion 5K ^ the quantity Q — 
{K)/{D — 1), cr — SK/\/D — 1, and ?/ is Gaussian noise 
with zero mean and unit variance. The factors involving 
D reflect the fact that the curvature-driven motion in 
the i*'' direction is, on the average, 'shared' by the D — 1 
orthogonal directions. 

Presuming further that the flow is ergodic and that (al- 
most) all orbits are chaotic, the quantities (K) and 5K 
can be calculated assuming a uniform sampling of the 
constant energy hypersurface, i.e., a microcanonical dis- 
tribution /i oc 5d{H — E). Alternatively, for sufficiently 
high dimensions one can proceed instead by assuming a 
canonical distribution which, for large _D, is much simpler 
computationally although it should yield nearly identical 
results. 

To complete the characterisation of the random pro- 
cess k{t) it remains to specify the autocorrelation func- 
tion T{t) or, at least, the autocorrelation time r, which 
governs how rapidly the curvature fluctuates along the 
orbit. This Casetti et al. provide using another geomet- 
ric argument. On the one hand, they identify a time 
scale 



Tl 



2^/TL- 



(4.9) 



corresponding to the typical time between successive 
conjugate points, i.e., points where the Jacobi field of 
geodesic deviation vanish. On the other, they identify a 
time scale 



T2 



(4.10) 



corresponding to the length scale on which the fluctua- 
tions become comparable to the average curvature. They 
then select as an appropriate autocorrelation time a value 
T satisfying 



4 



2 (rf^ 



so that 



2t : 



n 



2^/n{n + (T) + 7r(T 



(4.11) 



(4.12) 



Casetti et al. suggest further that T{t) might be well 
approximated by the osciUating function 



T(t) 



smut 



which yields an autocorrelation time 



1 

2w 



' osc 

1^' 



(4.13) 



(4.14) 



with Tosc the oscillation period. However, this is largely 
irrelevant for their analysis. Granted the assumption of 
additive Gaussian noise, the form of the color only enters 
into the final expression for x through the autocorrelation 
time T 

Given a knowledge of r and the first two moments, Eq. 
(4.8) can be solved analytically using a technique devel- 
oped by van Kampen [Q to yield an estimated largest 
Lyapunov exponent 



1 / 4^! 



where 



A = 



2crV 




1/3 



(4.15) 



(4.16) 



B. Applying this proposal to lower-dimensional 
Hamiltonian systems 

The preceding can be reformulated without recourse to 
differential geometry in a setting that makes the physi- 
cal content of the assumptions more transparent and, as 
such, makes it clearer which assumptions might prove 
suspect, especially for lower-dimensional systems. 

The basic perturbation equation (4.5) can be derived 
trivially from the original Hamilton equations 



dq' _ dH 
dt dpi 



and 



djh _ dH_ _ _dV 
dt dq^ dq'- 

associated with the Hamiltonian 

H = T + V =^5'^p,p,+Viq). 



(4.17) 



(4.18) 



(4.19) 



It is clear that the introduction of a small perturbation 
^ -(- and Pi ^ Pi + Ci leads to evolution equations 
of the form 



dt 



and — — 



dt 



dq^dq^ 



(4.20) 



but combining these last two equations leads immediately 
to Eq. (4.4). 

The crucial assumption underlying the entire Casetti 
et al. analysis is the assumption that, for the case of 
chaotic orbits, Eq. (4.4) can be modelled by a stochastic- 
oscillator equation. For the case of 'wildly chaotic' or- 
bits or orbit segments, which are far from periodic, this 
assumption would seem quite reasonable. However, in 
lower dimensional systems one encounters the possibil- 
ity of 'sticky' ||T^ orbit segments which, albeit charac- 
terised by positive short-time Lyapunov exponents, are 
'nearly regular' in visual appearance and have Fourier 
spectra with most of the power concentrated at or near 
a few special frequencies [|l6|. This is especially common 
for Z) = 2, where cantori ||l[| can serve as entropy bar- 
riers, confining a chaotic orbit near a regular island for 
surprisingly long times. To the extent that such orbit 
segments behave in a nearly regular fashion, the assump- 
tion of nearly random behavior is clearly suspect, and one 
might anticipate that a stochastic oscillator equation can- 
not prove completely satisfactory. Alternatively, to the 
extent that this 'sticky' behavior is rare, such an equa- 
tion might be expected to provide a reasonable starting 
point. 

The assumption of 'quasi-isotropy' can also be under- 
stood in very simple physical terms: Instead of consider- 
ing the Z?-dimensional Eq. (4.4), which involves the full 
second derivative matrix of V, it is assumed that, on 
the average, each direction in configuration space is sta- 
tistically identical, so that one can consider instead D 
identical one-dimensional equations. In this context, the 
only question concerns the proper identification of the 
quantity to play the role of the squared frequency k{t). 
The Casetti et al. prescription atates that the relevant 
information about stability is contained in the trace of 
the second derivative matrix, so that k{t) should be pro- 
portional to Ay = d'^V/dq'dq,. The factor of L> - 1 
entering into Eq. (4.8) reflect the fact that the perturba- 
tion driving the chaos is 'shared' amongst D — 1 dimen- 
sions. (Recall that, in a time-independent Hamiltonian 
system, there is always one direction of neutral stability 
corresponding to translation along the orbit from q'^{t) to 
q'{t + St).) 

This assumption of quasi-isotropy seems especially rea- 
sonable for large D where, on average, different directions 
of the configuration space should look much the same, 
but becomes more suspect in lower dimensions. In prin- 
ciple, one can relax this assumption by working with the 
full matrix equation. As a practical matter, however, this 
becomes quite cumbersome for D ^ 2. For this reason, 
most the following analysis will retain the assumption of 
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quasi-isotropy. What happens when this assumption is 
relaxed for two-dimensional systems is considered briefly 
in Section VI. It would in fact appear that, at least for 
D = 2, relaxing this assumption does not in general 
yield significant improvement in the estimated value of 
the largest Lyapunov exponent. 

For generic Hamiltonian systems with large D one an- 
ticipates that (almost) all the orbits on a constant-energy 
hypersurfacc arc chaotic. Granted the assumption of cr- 
godicity, it then follows that, over sufficiently long time 
scales, an orbit eventually samples a microcanonical dis- 
tribution. This implies that, when estimating a Lya- 
punov exponent x(-^)) ^ts defined in an asymptotic t ^ oo 
limit, one can assume that the statistical properties of 
the curvature experienced by an orbit are given correctly 
by a microcanonical distribution. By contrast, for lower 
dimensional systems, especially D = 2, one anticipates 
instead that a generic potential will admit a coexistence 
of large measures of both regular and chaotic orbits, so 
that the assumption of a microcanonical distribution is 
not justified. Rather, granted the assumption of ergod- 
icity, one would anticipate the existence of a different 
invariant distribution, corresponding to a uniform popu- 
lation of those portions of the constant energy hypersur- 
facc which are accessible to a chaotic orbit with specified 
initial condition. It is not clear how this distribution 
could be computed analytically. However, as described 
in Section V, numerical approximations to this invariant 
distribution can be generated straightforwardly through 
a time-series analysis of orbits evolved numerically. 

Even if a microcanical distribution is justified, the as- 
sumption of Gaussian fluctuations is problematic. For 
large D, this assumption can be motivated with a fair 
degree of rigor via a Central Limits Theorem argument, 
supposing that the distribution of values of AV^ involves 
a convolution of D nearly independent distributions for 
the separate components d^V / dq^dq^ (no sum over in- 
dices). For very small D, this is clearly not justified, and 
the minimum value of D for which the Gaussian approx- 
imation is justified must depend to a certain extent on 
the form of the individual distributions. It will be seen 
in Section V that, for the model systems (3.2) and (3.3), 
the convergence towards a Gaussian is quite efficient, and 
that the distribution A''[Ay] is reasonably well fit by a 
Gaussian even for D as small as 3 or 4. It will also be 
seen that, at least for distributions 7V[AF] that are not 
too skew, deviations from a Gaussian do not change the 
estimated value of the largest Lyapunov exponent all that 
much. 

The formula for the autocorrelation time r motivated 
by Casetti et al. is somewhat ad hoc in that, unlike the 
other crucial inputs fl and a, it cannot be derived di- 
rectly from a microcanonical distribution. However, the 
basic scaling implicit in t can again be inferred relatively 
simply. As will be seen below, fl and a are typically 
comparable in magnitude. They both reflect statistical 
properties of AV and, as such scale (within factors of 
order unity) asV/R"^, where V represents a typical value 



of potential and i? is a characteristic length scale, i.e., 
the size of the configuration space region probed by an 
orbit. Assuming 'virialisation,' i.e., that the mean po- 
tential and kinetic energies of the orbits are compara- 
ble in magnitude, it follows that V ^ v"^, where v de- 
notes a typical orbital speed. However, this implies that 
Q, a <^ v/R = f^, where to denotes a characteristic 
dynamical or crossing, time. Allowing for the fact that 
the characteristic scale on which V changes significantly 
is typically somewhat smaller than the size of region ac- 
cessible to the orbit leads to the obvious physical con- 
clusion that T should be comparable to, but somewhat 
smaller than, the time required for an orbit to travel from 
one side of the system to another. 

Implicit in the Casetti et al. analysis is the assump- 
tion that the stochastic process k{t) corresponds to state- 
independent, additive noise, so that, e.g., the autocorre- 
lation time r on which the curvature changes is indepen- 
dent of the value of the curvature. Strictly speaking, this 
assumption cannot be correct. If, e.g., an orbit is in a 
region where V is especially small, its kinetic energy, and 
hence its velocity, will be especially large, so that the or- 
bit will move very quickly to a different region where V , 
and hence in general AV, is very different. If, instead, 
the orbit is in a region where V is especially large, it 
will move more slowly so that Ay might be expected to 
change more slowly. The autocorrelation time r of (4.13) 
represents an average over a variety of orbits with very 
different values of V. One might anticipate that these 
differences will tend to wash out for large D, but there is 
no obvious reason why this should be true for smaller D. 

One final point: It is clear that, for small D, one cannot 
pass from a microcanonical to a canonical description. 
One must work directly with the microcanonical measure 
/i cx 5d{H — E). This, however, is not a major problem. 
For a D degree of freedom system, the microcanonical 
distribution corresponds to a configuration space density 

f iE - y)(^-2)/2 \iv <E- 
f{<f) cx <^ (4.21) 
[ i{V>E. 

but, given this formula for /, it is straightforward, at least 
numerically, to compute the distribution A/'[Ay] and/or 
any moments of the distribution. 

V. TESTING THE BASIC ASSUMPTIONS 
A. What was computed 

To test the validity of the basic assumptions requires 
a comparison of real orbital data with predictions made 
assuming a microcanonical distribution. The requisite 
orbital data were generated and analyzed as follows: 

For given choices of potential and total energy, a col- 
lection oi N = 1000 initial conditions were selected, and 
each of these was integrated into the future for a total 
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time T corresponding to between ^ 100 and 2000 char- 
acteristic crossing times tjj. The numerical integration 
simultaneously tracked the evolution of a small initial 
perturbation, periodically renormalised in the usual way 
so as to yield an estimate of the largest (short-time) 
Lyapunov exponent for the orbit segment. Configura- 
tion space data, recorded at fixed intervals 5t, were used 
to generate a time series {Kj{n5t)} for each of the seg- 
ments in the 1000 orbit ensemble which was deemed to 
be chaotic. 5t was typically so chosen that each seg- 
ment was sampled by 2560 points. Distinctions between 
regular and chaotic were implemented through the in- 
troduction of a threshhold value Xmm- if the computed 
X < Xmim the orbit segments were assumed to be reg- 
ular. Combining all the orbital data for all the chaotic 
orbits allowed the computation of the bulk moments {K) 
and 5K where, recall, K = /^V . Binning the combined 
data into 1000 bins yielded a numerical representation of 
the distribution N[K]. 

A discretized representation of the average autocorre- 
lation function T{t) was computed by selecting a rep- 
resentative ensemble of 5120 initial conditions, evolv- 
ing each of these into the future for an extended time 
T > 2048, so as to generate a set of well-mixed 'random' 
phase space points, identifying each of the Nc < N orbit 
segments that were chaotic, and, by extending the inte- 
grations for an additional time T' = 1024, constructing 

nn6t) = j^-^ J2 DK,{T)DK,{T + nSt). (5.1) 

Here DKj = Kj — (K), and the quantities (K) and {K"^) 
represent averages computed for all the chaotic orbital 
data for T < t < T'. Ideally one should compute the 
autocorrelation time t using the defining relation 

/•OC 

/ dtr{t) ^ {K^)t. (5.2) 

Given, however, that T is typically a rapidly oscillating 
function (period ^ to) with an envelope that damps very 
slowly, such a computation proved unreliable. A seem- 
ingly better measure of r or, at least, of how r scaled with 
energy E for fixed potential, was obtained by computing 
the period Tqsc associated with the oscillations. 

Predictions associated with a microcanonical distribu- 
tion were computed as follows: The microcanonical dis- 
tribution /i OC Sd{H — E) implies the configuration space 
probability density (4.22); but, given this /, it is straight- 
forward to compute the value of any configuration space 
function g{q). Numerical representations of the distribu- 
tion N[K] associated with a microcanonical distribution 
were computed by dividing the occupied configuration 
space into a collection of M hypercubes, (ii) deciding ran- 
domly whether or not to sample each hypercube, using 
a weighting (x {E — l/)(^~2)/2 ^ evaluated at a random 
point in the cube, (iii) in the event that the hypercube 
was to be sampled, locating a point in the cube at a 



randomly chosen location, and then (iv) binning the re- 
sulting collection of points into 1000 bins. 

Granted the assumption of a Gaussian distribution of 
curvatures, estimates of the Lyapunov exponent x can 
be, and were, computed using eq. (4.16), which does not 
require the assumption of a microcanonical population. 
When the assumption of a Gaussian distribution is re- 
laxed, an analytic solution is not possible in general, 
so that X was obtained instead from a numerical com- 
putation, with the random curvature generated initially 
by sampling A^[i4r], held constant for the autocorrelation 
time T, and then replaced by another, randomly chosen 
curvature pT| . 

B. What was found 

1. N[K] and its first two moments 

FIGURE ^ exhibits the energy-dependence of the 
quantities (K) and 6K for chaotic orbits in the dihedral 
potential with D = 2 and D = 3, computed both from 
time-series data (dashed lines) and assuming a micro- 
canonical distribution (solid lines). Overall, one observes 
excellent agreement between the numerical and analytic 
predictions, particularly for the first moment (K). The 
best overall agreement obtains for lower energies where, 
even for D = 2, the measure of regular orbits is compar- 
atively small and 'stickiness' seems comparatively unim- 
portant. 

For D = 2 at higher energies, say E > 1.0 or so, it 
appears that a third of the constant energy hypcrsur- 
face, or even more, corresponds to regular orbits, so that 
one is clearly not justified in assuming a microcanoni- 
cal distribution. However, it is evident that the predic- 
tions based on a microcanonical distribution remain quite 
good. That this should be the case is not really surpris- 
ing. Presuming that the regular islands are not concen- 
trated preferentially at regions where is especially 
large or small, it would seem reasonable to assume that, 
in a sufficiently coarse-grained sense, chaotic orbits still 
go 'all over' the energetically allowed regions of configu- 
ration space. To the extent, however, that this be true, 
one might expect moments approximating the moments 
appropriate for a microcanonical distribution which, for 
D = 2, implies [cf. Eq. (4.21)] a uniform configuration 
space density. FIGURE || exhibits analogous data for 
the FPU potential with D = 4 and D = 6, generated for 
parameter values a = 1.0 and b = 0.1. 

The thick solid curves in panels (a) - (d) of FIG. ^ ex- 
hibit distributions of curvatures, N[K], for the dihedral 
potential with D = 2 and D = Q generated assuming a 
microcanonical distribution. The corresponding curves 
in FIG exhibit analogous distributions for the FPU po- 
tential for Z? = 4 and D = 6. Panels (e) and (f) in FIG. | 
show the time series and microcanonical predictions for 
the dihedral potential in (from left to right) = 3, 4, and 
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5. It is evident that the microcanonical distributions for 
the dihedral potential with D = 2 are not even remotely 
Gaussian in shape. However, it is also apparent that, for 
all the other cases, the distribution is in fact reasonably 
well fit by a Gaussian, although N[K] typically has a 
slight skew and can manifest appreciable deviations for 
large \K - {K)\. 

The other curves in FIG. | (a) - (d) and in FIG. @ 
represent distributions TV [if] generated from time-series 
data. FIG. |^ and the first three panels of FIG. || display 
two numerical curves, one representing data for < t < 
1024 and the other for 2048 < t < 3072. FIG. | (d) also 
includes a third numerical curve, generated for 8192 < 
t < 9216. For the two energies exhibited in the D = 2 
dihedral potential, E — 1.0 and E — 6.0, there exist large 
measures of both regular and chaotic orbits and, for this 
reason, the time-series N[K] differs significantly from the 
microcanonical A'^[i^]. However, the ensembles of initial 
conditions used to generate the time-series distributions 
evolved towards an invariant (albeit non-microcanonical) 
distribution relatively quickly, so that the two numerical 
curves very nearly overlap. 

For the dihedral potential with D > 3 and for the FPU 
potential with D > A almost all the orbits appear to be 
chaotic, so that, assuming ergodicity, the microcanonical 
N[K] and a truly representative time-series N[K] should 
coincide up to statistical uncertainties. However, for the 
cases exhibited in FIGS. | (c) and (d) and FIG. @, the ini- 
tial ensembles only converged towards an invariant distri- 
bution on a comparatively long time scale, so that the two 
(or more) time-series curves differ appreciably from one 
another. In each case, the later time sampling(s) yielded 
distributions N[K] that more closely approximated the 
microcanonical A^[i^]. 

The preceding suggests that one can use the form of 
the distribution N[K] as a robust diagnostic in terms of 
which to probe the approach towards ergodicity. Ergod- 
icity per se is an assumption regarding the i — > oo limit 
and, even assuming ergodicity, there remains an obvious 
question: How long must one evolve some ensemble of 
initial conditions before its time-averaged density closely 
approximates the density associated with a constant pop- 
ulation of the accessible phase space regions? Comparing 
the distribution A^[_R'] associated with an evolving en- 
semble with the N[K] associated with a microcanonical 
distribution can provide a useful diagnostic for probing 
the extent to which the ensemble has evolved towards a 
time-independent invariant distribution. 

It is well known that different chaotic orbit segments in 
the same connected phase space region can exhibit vastly 
different short-time Lyapunov exponents, and that the 
values of these short-time exponents can correlate signifi- 
cantly with position. For example, chaotic segments near 
regular islands tend to be much less unstable than wildly 
chaotic segments located in the middle of the stochas- 
tic sea. One might, therefore, expect that orbit segments 
with especially large or small short-time exponents would 
be characterised by different curvatures. For potentials 



and energies where almost all the orbits are chaotic and 
'stickiness' is rare, this segregation should be minimal; 
but for potentials where there is a coexistence of large 
measures of both regular and chaotic orbits, and where 
'stickiness' is pronounced, this effect should be much 
more pronounced. 

As illustrated in FIG. ^, this intuition was corrobo- 
rated numerically. The top panel of FIG. || was gener- 
ated for E — —0.5 in the D = 2 dihedral potential, an 
energy where the regular regions are extremely small, so 
that a representative ensemble of 1000 initial conditions, 
integrated for a time T = 1024, yielded no regular orbits. 
The orbits generated from these initial conditions were 
divided into five quintiles, depending on the values of 
their short-time Lyapunov exponents, and the lower five 
curves in this panel exhibit individual subdistributions 
A^[i4'] computed for each quintile. The four quintiles cor- 
responding to the larger values of x yielded distributions 
that were nearly identical. The lowest quintile was again 
quite similar, but did manifest some noticeable differ- 
ences: This subdistribution, corresponding to the thick 
solid line, is distinctly underrepresented at very low val- 
ues of K and over-represented at very large K, and, un- 
like the other four quintiles, appears to be a slowly de- 
creasing function of K in the interval < K < 7.5. The 
sum of these five subdistributions (with a different nor- 
malisation from the subdistributions) corresponds to the 
slightly jagged upper curve, which is essentially identical, 
modulo statistical uncertainties, to the smoother curve 
computed for a microcanonical distribution. 

The lower panel of FIG || exhibits analogous data for 
E — 6.0, again in the D = 2 dihedral potential. In 
this case, a 1000 orbit ensemble was divided instead 
into a 'quintile' of 332 regular orbits and four 'quintiles' 
each comprised of 167 chaotic orbits, but the resulting 
subensembles were analyzed identically. The lower solid 
curve peaking at if ~ 13 represents the 332 regular or- 
bits, and the three nearly identical curves which have 
a local minimum at K ^ 13 correspond to the chaotic 
orbit segments with the largest short-time Lyapunov ex- 
ponents. The intermediate dashed curve corresponds to 
the chaotic orbits with the smallest short-time Lyapunov 
exponent, the majority of which could be reasonably clas- 
sified as 'sticky.' The total N[K] given as a sum of the 
four chaotic 'quintiles' is represented by the upper curve 
with a local minimum at if ~ 13. The upper curve cor- 
responding to a nearly fiat profile again corresponds to a 
microcanonical distribution. 



2. The autocorrelation time r 

As suggested by Casetti et ai, the autocorrleation 
function T{t) is in fact an oscillating function of time, 
but it tends to decay more slowly than with the 1/t en- 
velope implicit in Eq. (4.14). This slower decay is espe- 
cially evident for potentials and energies when 'stickiness' 
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is important, in which case a substantial 'memory' can 
persist for dozens of oscillations. This is, e.g., evident 
in FIG. ^ (a) - (d), which exhibit data for the dihedral 
potential for D = 2 and D = 6. The first two panels 
correspond to a very low energy E = —0.5, where, even 
for D = 2, almost all the orbits are chaotic. The second 
two panels correspond to a higher energy E = 6.0 where, 
for both D = 2 and D = 6, chaotic segments can be 
nearly periodic and have comparatively small short-time 
Lyapunov exponents. For D = 2, the case exhibited in 
panel (c), roughly one quarter of the chaotic orbits are 
noticeably 'sticky'; for D = 6, the case in panel (d), the 
fraction is reduced to about 5%. In either case, analy- 
sis of a sample that excludes 'sticky' segments yields an 
autocorrelation function that decays substantially more 
quickly. 

As noted by Casetti et ai, if the autocorrelation fimc- 
tion is in fact well approximated by Eq. (4.13) the time 
scale identified geometrically in Eq. (4.12) should co- 
incide with the time scale (4.14). This prediction was 
tested numerically and found typically to be satisfied to 
within factors of ~2, although some discrepancies were 
observed. As noted already, a direct determination of r 
using Eq. (5.2) proved unreliable. 

Perhaps the most striking point is that, at least when 
'stickiness' is comparatively unimportant, the Casetti et 
al. time scale t given by Eq. (4.13) and the time scale 
Tosc exhibit very similar scaling with energy E. This 
is illustrated in FIG. ^ (e) and (f), which exhibit both 
time scales as functions of E for the D = 2 and D = 6 
dihedral potential. In each case, the time scale Tosc is 
somewhat longer than the Casetti et al. time scale r. 
Significantly, though, for D ^ 2 the quantities Tqsc and r 
exhibit very different scalings at higher energies, precisely 
where 'stickiness' is most prominent. 

3. Sources of uncertainty 

Granted the assumption of a Gaussian distribution of 
curvatures, the predicted value of the largest Lyapunov 
exponent depends on only three quantities, namely (K), 
SK, and t; and as such, it is natural to ask how the pre- 
dicted value Xest varies if any of these inputs are changed. 

If one introduces a simulataneous scaling of both (K) 
and SK, i.e., (K) — s- a{K) and 6K — > aSK, with a of 
order unity, Xest ^ ct^^^Xest- If, alternatively, (fT) is held 
fixed but SK aSK , one infers, at least approximately, 
that Xest ct^^'^Xest- Finally, if (K) and SK are held 
fixed, but one allows for a scaling t ar, Xest ctXest- 

Given that the values of {K) and SK can be estimated 
quite well using simple dimensional arguments ~ recall 
that, even when there are relatively large measures of 
periodic orbits, a microcanonical population yields esti- 
mates in good agreement with the results of direct nu- 
merical computation -, it would seem that, with the as- 
sumption of quasi-isotropy and a Gaussian distribution 



of curvatures, the largest source of error should be in 
the determination of the autocorrelation time r. The ex- 
pression for r motivated by Casetti et al. is more ad hoc 
than the expression for (K) and SK; and dimensional ar- 
guments are hard pressed to yield estimates of r that are 
accurate to better than a factor of two. However, factors 
of two uncertainty in r translate directly into factors of 
two uncertainty in Xest- 

One can also investigate how the predicted Xest 
changes if one relaxes the assumption of a Gaussian dis- 
tribution, instead computing Xest by solving Eq. (4.9) nu- 
merically for the distribution [i^] generated either from 
a microcanonical distribution or from real orbital data. 
The resulting change in Xest will of course depend on the 
degree to which N[K] deviates from a Gaussian, larger 
deviations resulting in larger changes. Especially for two- 
dimensional systems, whre A'^[Ar] is far from Gaussian, 
allowing for the correct distribution can change Xest by 
factor of three, or even more. 

This is illustrated in FIG. ^, which exhibits several dif- 
ferent estimates of the largest Lyapunov exponents Xest 
for the D = 2 dihedral potential, most of which will be 
described in Section VI. In the present context, note sim- 
ply (i) the 'true' Xnum, generated by tracking a small ini- 
tial perturbation (solid line), (ii) the Casetti et al. Xest, 
based on the assumption of Gaussian fluctuations and an 
autocorrelation time given by Eq. (4.13) (short dashes), 
and (iii) an alternative Xest, again based on the 'quasi- 
isotropized' Eq. (4.9), but now allowing for a distribution 
Af[A'] generated from time-series data and an autocor- 
relation time (4.15) (long dashes). Both estimates are 
comparable in magnitude to Xnum, but both miss the 
nontrivial dip that arises near E — 0.0. 

VI. ESTIMATES OF LYAPUNOV EXPONENTS 
IN LOWER-DIMENSIONAL HAMILTONIAN 
SYSTEMS 

A. Estimates of the true Lyapunov exponent 

Overall, Eq. (4.15) first proposed by Casetti et al., 
modified to allow for moments computed assuming a 
microcanonical distribution, appears to give reasonable 
estimates of the largest Lyapunov exponent in lower- 
dimensional Hamiltonian systems, provided that an ap- 
preciable fraction of the phase space corresponds to 
chaotic orbits. In particular, as long as the true Lya- 
punov exponents are not very small {xnuml^ t^) and/or 
'stickiness' is not especially prominent, the estimated Xest 
typically agree with the numerical Xnum to within factors 
of two. In some cases, such as for the FPU model, the 
agreement between Xnum and Xest rapidly increases with 
increasing D; but in other cases, such as for the dihe- 
dral potential, this is not the case. This would suggest 
that the quasi-isotropy assumption, which should become 
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increasingly justified in higher dimensions, is not neces- 
sarily the principal source of error. 

FIGURE compares the numerical and estimated 
x{E) for the dihedral potential for = 2, 3, 4, and 
6. The estimated values were computed using Eq. (4.13), 
based on a Gaussian distribution, with moments gener- 
ated both from a time-series analysis (dashed lines) and 
assuming a microcanonical distribution (dotted lines). 
The numerical values are connected with a solid line. One 
observes significant differences in the shapes of the curves 
associated with the numerical and estimated values, but 
there is invariably an overall agreement to within a factor 
of two. The most striking discrepancies arise ioi D — 2, 
where the estimates completely miss the abrupt dip in 
Xnum for £^ 0. The fact that, for D = 2, the two esti- 
mated curves differ significantly at high energies reflects 
the fact that the constant energy hypersurface contains 
large regular islands, so that the invariant distribution is 
far from microcanonical. 

FIGURE |l2| exhibits the numerical and estimated 
X{E) for the FPU model for D = 4, 5, and 6, with 
the estimated values again computed assuming Gaussian 
distributions and moments generated from a time-series 
analysis. The data have been plotted on a log-log plot to 
facilitate comparison with FIG. 3 in Casetti et al. Here 
two points are immediately obvious: (1) The estimated 
and numerical curves are distinctly different, with Xest 
always larger than Xnum, but their curvatures are com- 
paratively similar. (2) The agreement between Xnum and 
Xest becomes progressively better for higher dimensions 
and for higher energies. For i? = 5 in = 4, where 
the numerical Xnum = 0.211 corresponds to a very long 
time t 45 3> in, Xest overestimates Xnum by nearly a 
factor of seven. For D = 6 and E = 5, Xest yields a value 
approximately 2.65 times too large; for E — 10240, its 
value is only 1.27 times too large. 

FIGURE ^ (a) exhibits the same data for the sixth 
order truncation of the Toda potential. As for the FPU 
model, Xest systematically overestimates the true Xnurm 
the largest errors arising at low energies, where Xnum 
is comparatively small, larger regular islands exist, and 
'stickiness' is especially important. 

B. Short-time Lyapunov exponents 

The computations described in the preceding subsec- 
tion indicate that, for a variety of lower-dimensional sys- 
tems, the Casetti et al. model of a stochastic-oscillator 
equation yields reasonable estimates of the largest Lya- 
punov exponent, x, as a function of energy E. How- 
ever, if the stochastic-oscillator picture is to be accepted 
as completely valid, one must also demand that it 'ex- 
plain' the varying degrees of chaos manifested by different 
chaotic orbit segments with the same energy, as probed 
by short time Lyapunov exponents. In particular, one 
might hope that, even if the estimated values Xest of the 



true Lyapunov exponent disagree significantly with the 
values Xnum computed numerically, the estimated and 
computed values of short-time Lyapunov exponents for 
different orbit segments with the same energy should be 
strongly correlated. For example, chaotic segments for 
which the true short-time Xnum is especially small should 
yield estimates Xest that are also especially small. 

That such correlations do in fact exist is illustrated in 
FIGS. ^ (b) and (c) and FIG. |l|, which exhibit results 
appropriate for, respectively, the truncated Toda and di- 
hedral potentials. Each of these FIGURES was gener- 
ated by (i) selecting a representative ensemble of 1000 
initial conditions, all with the same energy; (ii) evolving 
these into the future for a time T — 1024 while simulta- 
neously tracking the evolution of a small perturbation so 
as to generate Xnum', (hi) recording the values of K for 
each orbit at fixed intervals 5t = 0.4; (iv) analyzing each 
orbit to extract {K) and 5K; and (v) using these two 
moments along with Eq.(4.13) to generate an estimated 
Xest ■ The scatter plots provide unambiguous visual con- 
firmation that the values of Xnum and Xest are strongly 
related. 

This visual impression can be quantified by computing 
the Spearman rank correlation TZ between the values of 
Xest and Xnum in each ensemble, which satisfies 

6 ^ 

nXnum,Xest) = 1 - ^3 ^ ^ J2 ' (^'l) 

1=1 

Here TV = 1000 denotes the number of orbits in the en- 
semble and 5i denotes the difference in rank for the ith 
orbit when ordered in terms of Xnum and Xest . 7?, = 1 cor- 
responds to a perfect correlation; TZ — —1.0 corresponds 
to a complete anti-correlation. 

The data sets in FIG. |l^ (b) and (c), corresponding 
to = 30 and i? = 50 in the truncated Toda potential, 
both yield 7^ « 0.88. The data sets in FIG. |l| (a) - (e), 
corresponding to E — 1.0 in the dihedral potential, yield 
rank correlations ranging between a low of TZ ~ 0.85 for 
D = 3 and a high ofTZ ^ 0.95 for D = 2. The especially 
high rank correlation for D = 2 might seem surprising 
since the ensemble contains a large number of regular 
orbit segments, with very small Xnum- The reason TZ 
remains as large as it does is that, for orbit segments that 
are manifestly regular, so the Xnum eventually decays to 
zero, there is a correlation between the estimated value 
Xest and the rate at which Xnum tends to zero: for regular 
orbits where the short-time Xest is especially large, the 
convergence is especially slow, so that, at finite times, 
Xnum will also be especially large Q . 

Not surprisingly, the computed value of TZ for a given 
ensemble of initial conditions depends on the total in- 
tegration time. If the orbits be integrated for a suffi- 
ciently large T, their differences tend to 'wash out,' so 
that the observed range of values for Xnum and Xest both 
decrease. Eventually, the differences between different 
orbit segments become small and the pronounced corre- 
lation disappears. 
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It would seem visually from FIG. |TJ (a) - (e) that the 
numerical and estimated values of the short-time Lya- 
punov exponents deviate largely because of some overall 
scaling. Given that at least for D >3, the phase space is 
almost entirely chaotic, so that the distribution of curva- 
tures reflects a microcanonical distribution, the evidence 
(cf. FIG. ^) that this implies a nearly Gaussian distribu- 
tion, and the argument in the preceding subsection that 
quasi-isotropy is not necessarily the principal source of 
discrepancies, it would seem natural to conjecture that 
this reflects a misidentification of the proper autocorre- 
lation time T. Panel (f) in FIG. |lj shows what happens 
to the estimated value Xest for the D — 6 dihedral poten- 
tial if, for each orbit, r is reduced by a factor of w 0.75, 
the value required to ensure that, for the ensemble, the 
mean values of the estimated and numerical exponents 
coincide, i.e., (xnum) = {Xnum)- The net result is that, 
to a fair degree of approximation, the data points are 
aligned along Xest = Xnum- 



C. The special case D — 2 

It is natural to ask whether one can relax the assump- 
tion of quasi-isotropy, at least for D = 2 where it would 
seem most suspect. One way in which to do this would 
be to work instead with the Jacobi metric, which, for 
D = 2, leads to a single oscillator equation of the form 

i M 



dt^ W dt dt ^ ' 



(6.2) 



where W ^ E — V (q'^) denotes the kinetic energy and 



(6.3) 



Unfortunately, however, this equation is very difhcult to 
explore numerically since it yields near-divergences for 
W 0, which prove quite common for D = 2. 

Alternatively, within the setting discussed hitherto in 
this paper, there are two ways in which one might pro- 
ceed: 

1. Consider the full multi-dimensional Jacobi equation 
and view it as a matrix stochastic equation. This could 
at least provide the 'average' rate of exponential insta- 
bility in different configuration space directions. Quite 
generally, a small perturbation will satisfy 



(6.4) 



with V.,j = d^V/dq'dq 



the second derivative matrix. 
The objective then is to view each component of this ma- 
trix as an (approximately) independent stochastic vari- 
able, i.e., considering 



with Sriij a random variable. Given distributions A'^[V2;a:], 
7V[Vyy], and A^[T4y] — N\Vyx\, which can be computed 
from time-series data or assuming ergodicity, and some 
approximation to the autocorrelation functions Txx, ^yyi 
and Txy^ which can again be motivated cither from a time 
series or assuming ergodicity, this system is easy to solve 
numerically. |20[ 

2. At least for D = 2, it is easy to diagonalise the matrix 
equation (4.2) at any given instant so as to obtain the 
eigenvalues of the stability matrix. The corresponding 
eigenvectors will then satisfy equations of the form 



dt^ 



- i^±^± = 0, 



where the time-independent eigenvalues satisfy 



u;± 



Vyy)±\ [Vx 



yyj 



4K 



xy 



1/2 



(6.6) 



(6.7) 



sn 



(6.5) 



Viewing u!± — flo^±+6fl± as stochastic variables leads to 
a pair of decoupled oscillator equations which are easy to 
solve numerically. In general, one might anticipate that 
the smaller eigenvalue will correspond to the more rapid 
growth rate. 

These alternative were tested in detail for the D = 2 
dihedral potential. The principal results are summarised 
in FIG. which shows the numerical Xnum{E) (solid 
line) as well as estimated values Xest{E) generated in 
four different ways. The short and long dashed lines, 
discussed already in the preceding section, correspond to 
the isotropized equation (4.9), assuming a microcanoni- 
cal distribution (short dashes) or using inputs generated 
from orbital data (long dashes). The triple-dot dashed 
curve represents the values generated for the coupled os- 
cillator system and the dot dashed curve represents the 
values generated by solving Eq. (6.6) for a;_. All the 
estimated curves yield values Xest that agree with Xnum 
to within factor of two, but none seems especially better 
than the others. 

The hypothesis that chaotic behavior in lower- 
dimensional Hamiltonian systems can be modeled by a 
stochastic-oscillator equation would appear robust in the 
sense that different implementations all lead to predic- 
tions that yield at least rough agreement with numerical 
integrations. However, none of the alternatives consid- 
ered here would appear 'completely right.' It seems likely 
that, in very low-dimensional systems, the details mat- 
ter sufficiently that no universal prescription will yield a 
completely accurate prediction. 



VII. CONCLUSIONS AND DISCUSSION 

The results described in this paper strongly cor- 
roborate the intuition that chaotic motion in lower- 
dimensional Hamiltonian systems can be visualized as 
random, so that the average instability of chaotic orbits. 
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and hence the values of their largest Lyapunov expo- 
nents, can be derived from a harmonic oscillator equation 
with a randomly varying frequency. Modulo straightfor- 
ward modifications, technical rather than conceptual in 
nature, the approach introduced by Casetti et al. for 
higher-dimensional systems also works reasonably well 
for systems with dimensionality as low as I? = 2. In this 
sense, as suggested in the Introduction, one would ap- 
pear to have a clear new paradigm in terms of which to 
interpret the origins of chaos in lower-dimensional Hamil- 
tonian systems. 

The precise numerical values of x predicted using this 
analytic approach are somewhat less accurate in lower 
dimensions than they are for much larger D, but it re- 
mains true that, in general, this approach yields predic- 
tions that are correct to within factors of two. In princi- 
ple one might hope to do still better but, as a practical 
matter, this would seem difficult if not impossible. The 
'obvious' alternatives considered in Section VI C yield 
somewhat different predictions for the largest Lyapunov 
exponents. In some cases these predictions are somewhat 
better than those based on Eqs. (4.12), (4.15), and (4.16) 
but, overall, they seem neither appreciably better nor ap- 
preciably worse. This would suggest that the predictions 
based on these equations are comparatively robust, in 
the sense that minor modifications do not yield vastly 
different results. However, this might also suggest that 
there is no single, universal modification that one might 
introduce which would yield near-perfect agreement for 
all potentials and energies. In point of fact, this is hardly 
surprising. There is every reason to expect that details 
which should 'wash out' in higher-dimensional systems 
will remain important in lower-dimensional systems. A 
'thermodynamic' description of chaos should work best 
for systems with many degrees of freedom. 

In this context, two significant points should be 
stressed. 

(1) Even when the predicted values Xest{E) of the 
'true' Lyapunov exponent Xnum{E) are off by as much 
as factors of two, one observes strong correlations be- 
tween Xest{E) and Xnum{E) for different orbit segments 
with the same energy. Orbit segments for which the pre- 
dicted Xest is low tend to have small short-time exponents 
Xnum] and, similarly, segments for which Xest is high tend 
to have larger values of Xnum- The physics entering into 
Eqs. (4.12), (4.15), and (4.16) allows one to distinguish 
clearly between orbit segments that are 'wildly chaotic' 
in visual appearance and have especially large short-time 
exponents and 'sticky' segments that are nearly regular 
in appearance and have comparatively small short-time 
exponents. 

(2) The largest discrepancies between the predicted 
and numerically computed Lyapunov exponents occur 
invariably for those potentials and energies where large 
portions of the chaotic sea correspond to 'sticky' orbits 
manifesting nearly regular behavior, in which case Xest 
can be much larger than the 'true' Xnum- This is exactly 
what one would expect. If large portions of the stochas- 



tic sea are 'sticky,' an orbit will spend much of its time 
evolving in a nearly regular fashion, but it is clear that, 
while manifesting such near-regular behavior, its motion 
cannot be characterised as essentially random. Indeed, 
as discussed more carefully elsewhere 0, chaotic orbit 
segments for which {K) and 5K assume values close to 
those characteristic of regular orbits tend systematically 
to have very small short-time Lyapunov exponents. 

The principal difference between the approach devel- 
oped in this paper and the approach introduced by 
Casetti et al. is that the statistical properties of the 
mean curvature K are not derived assuming a canonical 
distribution. For a truly conservative system, a thermo- 
dynamic description must, strictly speaking, be based on 
a microcanonical distribution, and it is only for large D 
that one can approximate such a 'correct' description by 
a description based on a canonical distribution. More- 
over, for very low dimensions, especially D ~ 2, even a 
microcanonical distribution is clearly unjustified. A mi- 
crocanonical analysis is based on the assumption that 
the entire constant-energy hypersurface is chaotic, but 
for lower dimensions, non-integrable systems typically ex- 
hibit a coexistence of regular and chaotic regions, both 
with significant measure. A correct analysis must in- 
volve deriving the statistics of the curvature only in the 
chaotic phase-space regions, a task which seems difficult 
analytically but, given an assumption of ergodicity, is 
straightforward to implement via a time-series analysis. 

In part, this work concerning chaos and the phase 
mixing of chaotic orbits was motivated in the context 
of nonequilibrium systems comprised of a large num- 
ber of interacting particles. Examples of such sys- 
tems include self-gravitating systems, e.g., galaxies, and 
charged-particle beams governed by external focusing 
forces and internal Coulomb forces (space charge). For 
both these examples fast evolutionary time scales have 
profound consequences. For galaxies, they are an inte- 
gral part of the formation process For beams, they 
limit the degree to which an accelerator designer can pre- 
serve the beam quality, especially insofar as the evolution 
is irreversible p2| . 

As mentioned in Section II, one way to infer the fastest 
time scale is to consider the interaction of a single particle 
with the coarse-grained potential formed by all the other 
particles. The problem then reduces to one involving a 
low-dimensional Hamiltonian, and the obvious question 
is to what extent statistical arguments concerning the 
behavior of chaotic single-particle orbits can be invoked 
to simplify the analysis further. All the examples pre- 
sented herein suggest that time scales in low-dimensional 
Hamiltonian systems inferred via the statistical methods 
of Casetti et al. are typically valid within a factor of order 
two. They also suggest that uncertainties in the compu- 
tation of these time scales are principally associated with 
uncertainties in the autocorrelation time, and that these 
time scales are comparatively insensitive to the choice 
of the invariant measure that weights the statistical av- 
erages. More importantly, however, our examples rein- 
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force the idea advanced by Cerruti-Sola and Pettini p^ ] 
that rapid mixing originates from parametric instabihty 
due to positive-curvature fluctuations along the geodesic 
trajectories of the particles over the configuration-space 
manifold. Cerruti-Sola and Pettini conjectured that: 
"This mechanism is apparently the most relevant - and 
in many cases unique - source of chaoticity in physically 
meaningful Hamiltonians." The diverse set of examples 
presented here would seem to corroborate their conjec- 
ture. 

The statistical analysis, however, does carry some 
caveats. It generally 'predicts' fast exponential mixing in 
potentials that are known a priori to be integrable and 
thereby admit only regular orbits that can only mix secu- 
larly. Examples include spherically symmetric potentials 
for which Poisson's equation generates nonuniform den- 
sity distributions, and special triaxial potentials such as 
the Staeckel potentials |Q. Thus the analysis provides 
no information as to what criteria are necessary and suf- 
ficient to establish a preponderance of globally chaotic 
orbits; it merely hypothesizes their existence. Related to 
this deficiency is the failure of the analysis to account 
for 'sticky' chaotic orbit segments that, when present, 
will tend to slow down the mixing. Real systems may, 
however, mitigate these caveats. For example, external 
noise, even with very small amplitude, is known to add 
greatly to the efRciency of chaotic mixing by overcoming 
'stickiness' ||2^]. And the presence of localized irregular- 
ities that have been coarse-grained away may increase 
the chaoticity of the orbits. An important point, how- 
ever, is that graininess which manifests itself in binary 
particle interactions is not an example of such localized 
irregularities. Graininess establishes diffusion of an orbit 
away from the trajectory it would have in the smooth po- 
tential but, at least for nonchaotic orbits, this diffusion 
involves a secular, rather than exponential, divergence of 
trajectories ||2^ . 

Because it is based on the Eisenhart metric, the present 
treatment is restricted to stationary systems. However, 
with a Finsler metric, the geometric method can also in- 
corporate potentials that are explicitly time-dependent 
and/or velocity-dependent For example, recent 

work involving the Henon-Heiles potential resulted in 
a geometric measure of chaos over the associated Finsler 
manifold that was used for fast computation of the sys- 
tem's Poincare surface of section. If used with a coarse- 
grained potential, the Eisenhart metric includes no mech- 
anism for changing the particle energies. In principle, it 
can be included with a Finsler metric based on a time- 
dependent coarse-grained potential; however, the gener- 
alization also requires specifying a suitable invariant mea- 
sure for the nonequilibrium system [29|] . 

One final point should be noted. In writing this paper, 
the authors have deliberately adopted a tact somewhat 
complementary to that adopted by Casetti et al. Rather 
than focusing on the differential geometry of spaces ad- 
mitting an Eisenhart metric, the discussion has, to the 
extent possible, been couched in the language of conven- 



tional Hamiltonian mechanics. Such an approach serves 
to make the physical ideas underlying the general ap- 
proach more transparent physically and, it is hoped, will 
make the picture of chaotic motion as a random process 
comprehensible to a substantially larger audience. 
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FIG. 2. (a) Numerical (diamonds) and analytic (solid line) 
estimates of the largest Lyapunov exponent for chaotic or- 
bits evolved in the 3-dimensional galactic potential (2.1) ith 
A, &^ = 1, = H- A as a func- 
1.0. (b) The same for 
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tion of A for total particle energy E 
E = 0.6. 



FIG. 1. (a) Numerical (diamonds) and analytic (solid line) 
estimates of the largest Lyapunov exponent for chaotic or- 
bits evolved in the 3-dimensional galactic potential (2.1) with 
= 0.75, = 1.0, and = 1.25 as a function of black hole 
mass Mbh for total particle energy E = 1.0. (b) The same 
for E = 0.6. (c) E = 0.4. 




-2.5 -1.5 -0.5 0.5 

logio Mrh 



X 



0.3 
0.2 

0.1 

0.0 

0.3 
0.2 

0.1 
0.0 

0.3 

0.2 

0.1 
0.0 

-2.5 -1.5 -0.5 0.5 

logio 

FIG. 3. Same as Fig. 0, but with all orbits restricted to 
the 2-dimensional (x, y)-plane. 
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FIG. 4. (a) The mean curvature {K) for chaotic orbits 
in the D = 2 dihedral potential as a function of energy E, 
computed assuming a microcanonical distribution (solid line) 
and extracted directly from orbital data (dashed line), (b) 
The associated dispersion 5K. (c) {K) for D = 3. (d) 6K for 
D = 3. 
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FIG. 5. (a) The mean curvature (K) for chaotic orbits 
in the D — 4 FPU potential with a — 1.0 and 6 — 0.1 as a 
function of energy E, computed assuming a microcanonical 
distribution (solid line) and extracted directly from orbital 
data (dashed line), (b) The associated dispersion 5K. (c) 
(A') for D = 6. (d) 5K for D = 6. 
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FIG. 7. (a) The distribution of curvatures, N[K], for 
chaotic orbits with energy i5 = 20 in the = 4 FPU potential 
with a = 1.0 and b = 0.1, computed assuming a microcanoni- 
cal distribution (thick soHd Hne) and from orbital data for an 
ensemble evolved for times t = 1024 and t = 4196. (b) N[K] 
for D = 4 and E = 320. (c) N[K] for D = 6 and E = 20. (d) 
N[K] for D = 6 and E = 320. 



FIG. 6. (a) The distribution of curvatures, N[K], for 
chaotic orbits with energy E = 1.0 in the D = 2 dihe- 
dral potential, computed assuming a microcanonical distribu- 
tion (thick solid line) and from orbital data for an ensemble 
evolved for times t = 1024 and t = 3172. (b) N[K] for D = 2 
and E = 6.0. (c) N[K] for D = 6 and E = 1.0. (d) N[K] for 
D = 6 and E — 6.0. (e) A'^[J5'], as generated from time-series 
data for t = 3172, for E = 1.0 with £) = 3, £> = 4, and D = 5. 
(f) N[K] for E = 1.0 with D = 3, Z) = 4, and D = 5, now 
generated assuming a microcanonical distribution. 
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FIG. 8. (a) The distribution of curvatures, N[K], for a 
representative ensemble of 1000 orbits with energy E = —0.5 
in the D = 2 dihedral potential. The five lower curves rep- 
resent subdistributions, generated by dividing the ensemble 
into five quintiles based on the values of the short-time Lya- 
punov exponents for the orbits. The near-horizontal upper 
curve represents the distribution N[K] predicted by a mi- 
crocanonical distribution; the other, more jagged upper curve 
represents the distribution N[K] associated with the full 1000 
orbit ensemble, given by the sum of the five lower curves, (b) 
The same for E = 6.0. 
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FIG. 9. (a) The autocorrelation function T{t) for chaotic 
orbits in the D = 2 dihedral potential with E = —0.5. (b) 
The same for L) = 6 and E = —0.5. (c) The same for D = 2 
and E = 6.0. (c) The same for Z) = 6 and E = 6.0. (e) 
The Casettiet al. time scale r of Eq. (4.12) (solid line) and 
the time scale Tosc/'iTT of Eq. (4.14) (dashed line), for chaotic 
orbits in the D = 2 dihedral potential at different energies E. 
(f) The same for D = Q. 
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FIG. 10. Estimated values of the largest Lyapunov ex- 
ponent for the D = 2 dihedral potential as a function of en- 
ergy: Xnum generated from direct numerical integration (solid 
curve), the Casetti et al. value, generated assuming a Gaus- 
sian N[K] and autocorrelation time t given by Eq. (4.12) 
(short dashes); an estimate based on Eq. (4.8), but now using 
the N[K] generated from a time-series analysis and r given by 
Eq. (4.14) (long dashes); an estimate based on Eq. (6.6), using 
and r given by Eq. (4.14) (dot-dashes); and an estimate 
based on the coupled oscillator system, assuming Gaussian 
fluctuations and r given by Eq. (4.12). 
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FIG. 12. Estimated values of Lyapunov exponents for or- 
bits in the FPU model, generated from numerical integrations 
(solid lines) and estimated using Eq. (4.12). (a) For Z) = 4. 
(b) D = 5. (c) D = 6. 



FIG. 11. Estimated values of Lyapunov exponents for 
orbits in the dihedral potential, generated from numerical in- 
tegrations (solid lines) and estimated using Eq. (4.12). (a) 
For D = 2. {h) D = 3. (c) D = 4. (d) D = 6. 




FIG. 13. (a) Estimated values of Lyapunov exponents for 
orbits in the truncated Toda potential, generated from numer- 
ical integrations (solid lines) and estimated using Eq. (4.12) 
(b) Short-time lyapunov exponents computed using Eq. (4.12) 
(Xesi) and generated from numerical integrations (Xn«m) for 
E = 30.0. (c) The same for E = 50.0. 



FIG. 14. (a) Short-time lyapunov exponents computed 
using Eq. (4.12) (Xest) and generated from numerical integra- 
tions (xnum) for E = 1.0 in the D = 2 dihedral potential, (b) 
The same for D = 3. (c) D = 4. (d) D = 5. (e) D = 6. (f) 
Xnum for Z) = 6 contrasted with revised estimates Xest gen- 
erated by rescaling the time scale r of Eq. (4.12) by a factor 
of 0.75. 



